Astronomy & Astrophysics manuscript no. ms'v2 


©ESO 2012 


July 2, 2012 





A planetary system with gas giants and super-Earths around the 

nearby M dwarf GJ 676A 

Optimizing data analysis tecliniques for tlie detection of multi-planetary systems 

Guillem Anglada-Escude^ , Mikko Tuomi^'^ 
' Universitat Gottingen, Institut fiir Astrophysik, Friedrich-Hund-Platz 1, 37077 Gottingen, Germany 

- University of Hertfordshire, Centre for Astrophysics Research, Science and Technology Research Institute, College Lane, ALIO 
(psj". 9 AB, Hatfield, UK 

^ University of Turku, Tuorla Observatory, Department of Physics and Astronomy, Vaisalantie 20, Fl-21500, Piikkio, Finland 

o 

(N 

C '. ABSTRACT 

Context. Several M dwarfs are targets of systematical monitoring in searches for Doppler signals caused by low mass exoplanet 
\ companions. As a result, an emergent population of high multiplicity planetary systems around low mass stars starts to show up as 

well. 

Aims. We optimize classic data analysis methods and develop new ones to enhance the sensitivity towards lower amplitude planets 
in high multiplicity systems. We apply this methods on the public HARPS observations of GJ 676A, a nearby and relatively quiet M 
I dwarf with one reported gas giant companion. 

, Methods. We rederive Doppler measurements from public HARPS spectra using the recently developed template matching method 

(HARPS-TERRA software). We use refined versions of periodograms to assess the presence of additional low mass companions. We 
also analyse the same dataset using Bayesian statistics tools and compare the performance of both approaches. 
Results. We confirm the already reported massive gas giant candidate and the presence of a long period trend in the Doppler 
Q , measurements. In addition to that, we find very secure evidence in favour of two new candidates in close-in orbits and masses in the 

< I ■ super-Earth mass regime. Also, the increased time-span of the observations allows the detection of curvature in the long period trend 

^ ■ suggesting the presence of a massive outer companion whose nature is still unclear. 

^ I Conclusions. Even though the increased sensitivity of our new periodogram tools, we find that Bayesian methods are significantly 

more sensitive and reliable in the early detection of candidate signals but more works is needed to properly quantify their robustness 
against false positives. While hardware development is important in increasing the Doppler precision, development of data analysis 
techniques can help revealing new results from existing data sets with significant less resources. These new planetary system holds 

>: 

OO ■ exoplanetary system with a general architecture similar to our Solar System. GJ 676A can be happily added to the family of high 

multiplicity planetary systems around M dwarfs. 
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■ 1. Introduction The list of planets detected by HARPS is long an varied as 

— ' can be seen in the 35 papers of the series The HARPS search 

Doppler spectroscopy IS currently the most effective method m ^^^^^^^^ extra-solar planets. Instead of citing all of them. 



X 



detecting planet candidates orbiting neai'by stai's. The current ^^^^^ the intereste d reade r to the lates t HARP S results pre- 
precision enables the detection of planets of a few Earth masses ^^^^^^ iPepe et all (l20Tlh : lMavoretal.1 (l20Tlh : iBonfils et al.1 



5, • i£^se-in orbits, especially around low mass stars (|Mayor et alj ^^^^ demonsti-ated a radial velocity (RV) stabil- 

^ ■ 120091). Two methods are currendy used to get precision Doppler j ^j^^ ^^^^^ ^ ^ ^-1 jij^e-scales of several years. Since 

■ ■ measurements in the visible part of the spectrum : the gas j^^^^ 2OII, reduced data products derived from the HARPS 

cell and the stabilized spectrograph approach The Iodine cell j^^^^ Reduction Software (DRS) are publicly available through 

method consists on inserting a cell containing Iodine gas in the ^ dedicated webpage in ESO Q All data used in this work has 

beam of the telescope thus providing a very accurate method to ^^^^ obtained from there 

solve for the wavelength solution, instrumental profile variabil- r,,, • • ■ j ^ 1 • 1 1 • • 1 

, ^, „ , 1° • . 11 . 1^T^T^ r~n The increasing demand for higher Doppler precission has 

ity and the Doppler changes in the stellar spectrum OButler et alj , . .J' . • , ■ , , ■ 



, J u • 1- J T. Tj- Z — ■ motivated a Significant investment in hardware development and 

199a). The second approach is based on building a mechani- , ^ , ... , , ■ .1 

— TT"^ . 1 1 £i_ r J ^ 1 1-1 J • 1 • • a number of new stabilized spectrographs are currently un- 

cally stable fiber-ted spec trograph calibrated with an emission , . / ,^rii . 1 ^^v1^^ • 1 u 

lamj Raranne et al Il l99fi h HARPS is the best examnle of a sta COnstructlon(e.g., Wllken et al. 2012). It IS known, how- 

laiiip upgranne ei ai. lyy u*). n/vivro is Liie uesL exaiiinie 01 a SLa- , , , , . 1 1 ,1 t» a T-»T-.r~i T^T-»r~i 

u-i- ^ 1 J • . 11 J . .u o ,c ever, that the method employed by the HARPS-DRS to ex- 

biuzed spectrograph in operation and It IS installed at the 3.6m ' , , ■' , ■ ^ • ^• 1 

„, . T T-n^^ 1 -1 Jr. i ii^rvrvok tract RV mcasurcmcnts (Cross-correlation tunction) IS subopti- 

Telescope in La Silla-ESO observatory/Chile (IPeoe et aIJI2003l) . , . ... , , ■ ■ , 

^ ■' — ^ mal in exploiting the Doppler information in the stellar spec- 
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tram (IOueloziri995t iPepe et aT]l2002h . so developments in the 
data analysis techniques are also required to reach photon 
noise limited precision. We have recently developed a least- 
squares template matching meth o d (HARPS-TERRA software , 
lAnglada-Escude & Butlei I20TI lAnglada-Escude et all I20T2I) 
that is able to derive precise RV measurements from HARPS 
spectra obtaining a substantial increment of precision on K and 
M dwarfs. In (Anglada-Escude & Butler 2012), 34 observations 
on the MO dwarf GJ 676A were used to illustrate the improve- 
ment in precision of the HARPS-TERRA measurements com- 
pared to those used i n the discovery paper of the massive planet 
candidate GJ 676b ( Forvei lle et al. I I2OIII) . Additional observa- 
tions on this star have been recently released through the ESO 
archive, so we applied HARPS-TERRA to extract new RV mea- 
surements of the full set. In a quicklook analysis using clas- 
sic periodogram methods, we found tentative evidence for addi- 
tional planets in the system. However, these preliminary detec- 
tions did not provide convincing results. Recent developments 
in the Bayesian analysis methods of Doppler data (e.g., iTuomil 
I2OI2I) indicate that correlation between parameters seriously af- 
fect the sensitivity of periodogram based methods in the de- 
tection of additional low amplitude signals. Moreover, careful 
Bayesian anal ysises provide increased sensitivity to lower am- 
plitude signals (lTuomill20l"2l) and seem to be less prone to false 
positives than methods based on sequential periodogram analy- 
ses of the residuals only(Tuomi 2011). 

In this work, we develop and test data analysis methods for 
optimal detection of low-mass companions in multi-planetary 
systems and apply them to the HARPS -TERRA measurements 
of GJ 676A. In Section|2l we describe a new periodogram-based 
approach (recursive periodogram) and review the Bayesian anal- 
ysis tools also developed to deal with multi-Keplerian fits. 
Section[3]reviews the stellar properties of GJ 676A, describes the 
observations, discusses periodicities detected in activity indices 
and describes the previously det ected candidates (one massive 
gas giant and a long period trend iForveille et aDl201 ih . Section 
|4] analyses the RVs of GJ 676A using the aforementioned tools. 
Both methods (recursive periodograms and Bayesian analyses) 
agree in the detection of two additional sub-Neptune/super-Earth 
mass candidates in close-in orbits. We also use the opportunity to 
test the sensitivity of both detection methods by applying them to 
a subset of observations (first 50 epochs). We found that, while 
the recursive periodogram approach could only spot one of the 
signals at low confidence, a Bayesian analysis could already re- 
cover (or strongly suggest) the same 4 candidates much earlier. 
Finally, in Section|5]we place the unique features of the planetary 
system around GJ 676A in the context of the currently known 
population of exoplanets and provide some concluding remarks. 



2. Data analysis methods 

2.1. Recursive periodograms 

Classic least- s quares periodog rams and derived methods 
(IScargld IT982L |C ummind l2004l) consist of adjusting a sine- 
wave (equivalent to a circular orbit) to a list of test periods 
and plot these periods against some measure of merit. When 
k-periodic signals are detected in the data, the corresponding 
Keplerian model is subtracted from the data and a least-squares 
periodogram is typically applied on the residuals to assess if 
there is a k+lth periodicity left. As noted by several authors 
dAnglada-Escude et al. 2010; Lovis et al. 2011b; Tuomi 2012), 
non-obvious correlations between parameters are likely to de- 
crease the significance of (yet undetected) low amplitude signals. 



This is, as the number of the Keplerian signals in a model in- 
creases, the aliases of previously detected signals and other non- 
trivial correlations seriously affect the distribution of the resid- 
uals and, unless the new signal is very obvious, a periodogram 
on the residuals will not properly identify (even completely con- 
fuse) the next most likely periodicity left in the data. 

To account for parameter correlation at the period search 
level, we have developed a generalized version of the classic 
least-squares periodogram optimized for multiplanet detection 
that we call recursive periodogram. Instead of adjusting sine- 
waves to the residuals only, a recursive periodogram consists on 
adjusting all the parameters of the already detected signals to- 
gether with the signal under investigation. Even in the presence 
of correlations, candidate periods will show prominently as long 
as the new solution provides a net improvement to the previ- 
ous global solution. In our approach and by analogy to previous 
least-squares periodograms, a circular orbit (sinusoid) is always 
assumed for the proposed new periodicity. When no previous 
planets are detected, this is equivale nt to the GeneraUzed least - 
squares periodogram discussed by ' Zechmeister et alJ (l2009l) . 
and is a natural generalization of the methods discussed by 
iCummind (|2004|) to multi-Keplerian solutions. The graphical 
representation of the periodogram is then obtained by plotting 
the obtained period for the new planet (X-axis) versus the F- 
ratio statistic obtained from the fit (Y-axis). The highest peak in 
this representation indicates, in a leasts-squares sense, the most 
likely periodic signal present in the data. 

As any other classic least-squares periodogram method, one 
has to assess if adding a new sig nal is justified gi ven the improve- 
ment of the fit. As proposed by Gumming (2004), we use the F- 
ratio statistic to quantify the improvement of the fit of the new 
model (k+ 1 planets) compared to the null hypothesis (k planets). 
The F-ratio as a function of the test period is defined as 



F(P) 



xlliNobs-Nk^i) 



(1) 



where A^^ is the number of free parameters in the model with k- 
planets and Nk+\ is the number of free parameters in the model 
including one more candidate in a circular orbit at period P. For 
a circular orbit, the number of additional parameters Nk+\ - Nk 
is 2 (amplitude and phase of a sinusoid). Under the assumption 
of large Nobs, statistical independence of the observations, and 
Gaussian errors; F{P) would follow a Fisher F-distribution with 
Nk+\ - Nk and N„bs — Nk+\ degrees of freedom. The cumulative 
distribution (integral from to the obtained F-ratio) is then used 
as the confidence level c at each P (also called single frequency 
confidence level). Because the period is a strongly non-linear pa- 
rameter, each peak in a periodogram must be treated as an inde- 
pendent experiment (so-called independent frequencies). Given 
a dataset, the number of independent frequencies M can be ap- 
proximated by P,„inAT, where A T is the time baseline of the 
observations and P,„,„ is the shortest period (largest frequency) 
under consideration. Given M, the analytic false alarm probabil- 
ity is finally derived as FAP= \ - c'^ . 

Since the fully Keplerian problem is very non-linear, sev- 
eral iterations at each test period are necessary to ensure con- 
vergence of the solution. A typical recursive periodogram con- 
sists in testing several thousands of such solutions and, there- 
fore, it is a time consuming task. As a result, special care has 
to be taken in using a robust and numerically efficient model to 
predict the observables required by the solving algorithms. We 
found that a slight modification of the parameterization given by 
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(IWright & Howardll2009l) provided the best match to our needs. 
The only change we applied was using the initial mean anomaly 
Mo instead of the time of periastron To as a free parameter. These 
two quantities are related by IttTq/P = -Mo. From this expre- 
sion one can see that the replacement of rQ by Mp eliminates 
the n on-linear coupling between Tq and P. dWright & Howard 
12009 ) also provides the partial derivatives of the observables (ra- 
dial velocity) in a numerically efficient representation. The par- 
tial derivative of the RV with respect to Mq (instead of Tq) is 
trivially obtained as minus the partial derivative of the radial ve- 
locity with respect to the mean anomaly E {dv/dM^ - -dv/dE). 
Beyond this change, the adopted model is identical to the one 
given in (Wright & Howard 2009) so, for the sake of brevity, we 
do not provide the full description here. To accelerate conver- 
gence, at each test period, we first solve for linear parameters 
only. Next, we use the Levenverg-Marquadt(Levenberg 1944) 
method to smoothly approach the minimum and, finally, 
a few interations of a straight non-linear least-squares solver 
(iPress et al.lll992h are used to quickly converge to the final so- 
lution. Note that, although fitting for k-planets at each test pe- 
riod would seem a very time-consuming effort, we are implicitly 
assuming that the solver is already close to the local mini- 
mum, so a relatively low number of non-linear iterations (be- 
tween 20 and 50) are typically enough to reach the closest local 
minimum. Since all the orbits are re-adjusted and even though 
the method still suffers from some of the typical pitfalls of se- 
quential Keplerian fitting (e.g., one can still get stucked on local 
minima), the solution at each test period always has a higher sig- 
nificance than a periodogram on the residuals, especially in the 
presence of correlations between parameters. 

It is known that the assumptions required by the F-ratio tests 
might not be stricly satisfied by RV observations. Therefore, an 
empirical scheme is always desirable to better asses the FAP of a 
new detection. Since the recursive periodogram is a straight gen- 
eralization of least-squares periodograms, we adopt the brute- 
force Monte Carlo method proposed bv . Gumming (2004.) to ob- 
tain empirical estimates of the FAPs. This is, we compute re- 
cursive periodograms on synthetic datasets and count how many 
times spurious peaks with higher power than the signal under 
investigation are obtained by an unfortunate arrangement of the 
noise. Each Monte Carlo trial consists on : 1) taking the resid- 
uals to the model with k-planets and randomly permutate them 
over the same observing epochs, 2) adding back the signal of the 
model with k-planets, 3) re-adjusting the solution with k-planets 
(new null hypothesis), 4) computing the recursive periodogram 
on this new synthetic dataset and, 5) recording the highest F-ratio 
in a file. The FAP will be the number of times we get an F-ratio 
higher than the original one divided by the number of trials. 

A recursive periodogram can take a few tens of minutes de- 
pending on the number of datapoints and number of planets in 
the model. While this is not a serious issue while exploring one 
dataset, it becomes a problem when FAPs have to be empirically 
computed on many thousands of Monte Carlo trials. As a gen- 
eral rule, we accept new candidates if they show an empirical 
FAP lower than 1%. While this threshold is arbitrary, it guaran- 
tees that even if some of the proposed candidates are false pos- 
itives, we will not seriously contaminate the current sample of 
~700 RV detected exoplanets with spurious detections. As a first 
saving measure and given that analytic FAPs are known to be 
over-optimistic, we will only compute empirical FAPs if the an- 
alytical FAP prediction is already lower than 1%. The chance of 
obtaining a false alarm in a trial is a Poisson process and, there- 
fore, the uncertainty in the empirical FAP is ~ V-^FAPs/Mriais- 
Our aim is to guarantee than the empirical FAP is < 1 % at a 4-cr 



level, so we designed the following strategy to minimize the the 
number of Monte Carlo trials. This is, we first run 1000 trials. 
If no false-alarms are detected, the candidate is accepted and the 
analytic FAP is used to provide an estimate of the real one. If 
the number of trials generating false alarms is between 1 and 20 
(estimated FAP~0. 1-3%), we extend the number of trials to 10"*. 
If the updated FAP is lower than 0.5%, we stop the simulations 
and accept the candidate. If the empirical FAP is still between 
0.5 and 1.5%, 5 10"* trials are obtained and the derived FAP is 
used to decide if a candidate is accepted. While the computation 
time for 1000 trials in a single processor is prohibitively high, 
the computation of many recursive periodograms can be easily 
parallelized in modern multi-processor desktop computers. For 
the GJ 676 A dataset and a (3H-l)-planet model, 10^ trials would 
take 2.3 days on one 2.0 GHz CPU. The same computation on 40 
logical CPUs takes 1 .4 hours allowing to obtain empirical FAP 
runs with 10'*-10^ trials in less than a week. 

2.2. Bayesian analysis methods 

As in e.g. lTuon ii|(|2012|), the Bayesian analyses of the RVs of GJ 
676A were conducted using samplings of the posterior probabil- 
ity densities, estimations of Bayesian evidences and the corre- 
sponding model probabilities based on these samples. 

We sampled the pos terior densities u sing the adap- 
tive Metr opolis algorith m (iHaario et al.l 1200 ll) . also described 
broadly in lTuomil (1201 ih . Because it converges reliably and rela- 
tively rapidly to the posterior density in most situations, we per- 
formed several samplings of the parameter space of each model. 
Different samplings were started with different initial states to 
make sure that the global probability maximum of the param- 
eter space was found for each model. If they all converged to 
the same solution, we could confidently conclude that the corre- 
sponding maximum was indeed the global one. This check was 
performed because it is possible that the Markov chains get stuck 
in a local maximum if it is sufficiently high and the initial state 
happens to be close to it. As a result, we could then reliably 
estimate the parameters using the maximum a posteriori (MAP) 
estimate s and Bayesian credibility sets (BCSs) as uncertainty es- 
timates dTuomi & Kotirantall2009t) . 

The Bayesian evidences of each model were calculated 
using the one block Metropolis-Hastings (OBMH) estimate 
(IChib & Jeliaz kov 2001). It requires a statistically representative 
sample from the posterior, available due to posterior samplings, 
and can be used to assess the evidences and the correspond- 
ing model probabilities with relatively little computational effort 
when determining the number of Keplerian signals in a RV data 
set fa voured by the data (e.g. lTuomill201 ll 120121 : iTuomi et al.l 
120111). 

Using the OBMH estimates, we determined the probabili- 
ties of the models with differing numbers of Keplerian signals. 
However, we did not blindly choose the model with the great- 
est posterior probability and add it to the solution unless three 
detection criteria were also satisfied. We required that (1) the 
posterior probability of a model with k+l Keplerian signals was 
at least 150 times greater than that of a model with k signals 
(Kass &Raftery 1995; Tuomi 2011, 2012; Tuomi et al. 201 1|); 
(2) the RV amphtudes o f every signal were significantly greater 
than zero (lTuomill2012h : (3) and that the periods of each signal 
were well constrained from above and below because if this was 
not the case, we could not tell whether the corresponding signals 
were indeed of Keplerian natu re and periodic ones. These detec- 
tion criteria have been used in lTuomil ( 1201 2|) and they appear to 
provide reliable results in terms of the most trustworthy number 
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Table 1. Basic parameters of GJ 676A 



Parameter 


Valtie 


Reference 


R.A. 


17 30 11.203 


(a) 


Dec. 


-51 38 13.104 


(a) 




-260.02 ± 1.34 


(a) 


//Dec [mas yr"'] 


-184.29 ±0.82 


(a) 


Parallax [mas] 


60.79 ± 1.62 


(a) 


V 


9.585 ± 0.01 


(b) 


K* 


5.825 ± 0.03 


(c) 


Sp. type"^ 


MOV 


(b) 


Mass [MaY 


0.71 ± 0.04 


(d) 


Fe/H 


+0.23 ± 0.10 


(e) 



Notes , (a) HIPPARCOS catalog, (van Leeuwen 2007) (b) dKoen et alj 
l201Ch (c^ 2MASS catalog, (iSkrutskie et al.. .2006) (d) Using 
jPelfosse et al.ll2000l) (e) Using dJohnson & A'pDsll2009h 

of signals in a RV data set. We claim a detection of a Keplerian 
signal in the data if the Markov chains of several samplings con- 
verge to a solution that satisfies the criteria 1-3 above. 

The prior probability densities were chosen to have the same 
quantitative forms as in (Tuomi 2012), in which e.g. the param- 
eter space of the RV amplitude was limited to [0, 20] ms"'. 
However, because the RV data contains the obvious Keplerian 
signal of a massive ca ndidate and a long-period trend reported 
dForveille et al.l 1201 ll) with amplitudes clearly larger than 20 
m , the first two signals were allowed to explore a wider 
range o f semi-major amplitudes, i.e., [0, 200] ms Also, fol- 
lowing (lTuomill2012l ). we did not set the prior probabilities of 
different models equal but set them such that for models Mk 
and Aik+i, it holds that the prior probabilities satisfy P{Aik) - 
2P(Mk+i) for all values of k. 

3. Observations, stellar activity and previous work 

GJ 676 is a common proper motion pair of M dwar fs. The pri- 
mary (GJ 676A) has been classified as an MOV star ('Koen et al.l 
I2OIOI) . Using the empirical relat ions of [P elfosse et al. (2000), 
the 2MASS JHK photometry ( Skrutskie et al.i 1200 6) and its 
trigonometric parallax (.van Leeuwen 2007r we derive a mass 
of 0.71 Mq for GJ 676A. The star does not show strong evi- 
dence of activity nor evidence of youthness and, theref ore, it is 
a go od candidate for high precision RV studies(Forveil le et al.l 
I2OII ). The basic parameters of GJ 676A are given in Table [T] 
The fainter member of the pair (GJ 676B) has been classified 
as an M3V and is currently separated 50 " from A. From its 
HIPPARCOS parallax, this coiTesponds to a minimum separa- 
tion of 800 AU and an orbital period longer than 20 000 years. 
At this separation, the maximum acceleration of GJ 676A caused 
by GJ 676B on our line of sight is about 0.05 m s ' yr 

New radial velocity mea surements were obtained using 
the HARPS-TERRA software (Anglada-Escude & Butleij l2012l) 
from HARPS spectra recently made public through the ESO 
archive. The spectra are provided extracted and wavelength cal- 
ibrated by the HARPS-Data Reduction Sofwtare (DRS). Each 
HARPS spectrum consist of 72 echelle appertures covering the 
visible spectrum between 3800 and 6800 A. The average spec- 
tral resolution is A/SA - 1 10000 and each echelle apperture con- 
sists on 4096 extracted elements (or pixels). The set of public 75 
spectra have been obtained by several programs over the years 
and typical exposure times vary between 300 to 900 seconds. 
The mean signal-to-noise ratio (SNR) at 6000 Ais 60 and, in a 
few cases, it can be as low as 22. Doppler measurements derived 



with HARPS-TERRA are diff'erential against a very high SNR 
template spectrum generated by coadding all the observations. 
The secular acceleration eff'ect (Zechmeister et al. 2009) was 
subtracted from the RVs using the HIPPARCOS dvan Leeuw-enl 
l2007i) proper motion and parallax of the star. 

Two sub-stellar co mpanions were reported on the system by 
iForveille et al.l (|201 ih . The most prominent one is a massive gas 
giant candidate with a period of ~ 1060 days and a semi-major 
amplitude of ~ 120 m s ' . Strong eviden ce for a second, very 
long period candidate was also reported in iForveille et al.l(l201 ih 
be cause of a strong tren d detected in the residual to the 1 planet 
fit. IForveille et al.l (1201 1) akeady noted that the magnitude of 
such trend (~ 8 m s"' ) was too high to be explained by the grav- 
itational pull of GJ 676B (~ 0.05 m s~' ). Even after subtract- 
ing a model with one planet and a trend, ^ Forveille et al.l (1201 Ih 
also noted that the RMS of the residuals was significantly higher 
(~3.6 m s"') than the reported uncertainties (1 to 1.5 m s"' ) 
which was suggestive of pot ential new candidates Reanalysis 
of the 40 spectra available to lAnglada-Escude & Butleii (120 12h 
confirmed that, even with the increased precision derived using 
HARPS-TERRA (RMS 3.2 m s '), the star did show a signifi- 
cant excess of RV variability. 

In addition to the RV measurements, HARPS-TERRA also 
measures the Call H-hK activity in dex (S-index in the Mount 
Wilson svstem iBahunas et a l. 1995) and collects the measure- 
ments provided by the HARPS-DRS on two of the cross- 
correlation function (CCF) parameters that are also sensitive 
to stellar activity : bisector span (or BIS), full-width at half- 
maximum of the CCF (or FWHM). The S-index is directly mea- 
sured by HARPS -TERR A on the blaze-cor rected spectra using 
the definitions given by iLovis et aP (1201 la) and is an indirect 
measurement of the chromospheric emission and the average in- 
tensity of the stellar magnetic field. Because the strength of the 
magnetic field aff'ects the efficiency of convection, some spu- 
rious RV signals should correlat with the variability in the S- 
index. BIS is a measure of the asymmetry of the average spec- 
tral line and should coiTelate with the RV if the ob served ofi'sets 
are caused by spots or plages rotating with the star (Ouelo z et al.l 
2001). The FWHM is a measure of the width of the mean spec- 
tral line and its variability is usually associated with changes in 
the convective patterns on the surface of a star that might also in- 
duce spurious RV offsets. Since the connection betw een activity 
and RV jitter on M dwarfs is not well understood yet (iLovis et al.l 
I2OI lal) . we restrict our analyses to evaluate if any of the indices 
has periodicities similar to the detected RV candidates. 

No strong periodicities were detected on the BIS and the 
FWHM. However, the S-index did show a very strong signal 
with a period of 910 days (FAP~ 0.01%). Although the sig- 
nal has a similar period as the G J 676Ab candi d ate, su ch co- 
incidence was not mentioned in IForveille et al.l (l201 ll) . Even 
if the periods are similar, two strong arguments favor the 
Keplerian interpretation of the candidate. First, the signal in the 
S-index signal is not in pha se with the Doppler one. Second, 
iGomes da Silva et al.l (1201 2b has shown that RV offsets coiTe- 
lated with the variability of the S-index (or similar spectroscopic 
indices) are at the level of a few m 5"' level at most (GJ 676Ab's 
RV semi-amplitude is 120 m s"' ). Therefore, we find it likely 
that the similarity in the periods is purely coincidental. Evidence 
for an activity cycle of ~ 1000 d ays is also supported by the 
analysis of the Nal-index done by iGomes da Silva et al.l (1201 2|) 
using almost the same set of observations. 

In a preliminary analysis of the new 75 RVs, periodograms 
on the residuals to the two Keplerian solution (gas giant + trend) 
showed several tentative high peaks at 36, 59 and 3.6 days. 
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Fig. 1. Detection periodograms from most significant signal to 
less significant one (top to bottom). Black lines are least-squares 
periodograms computed on the residuals to the k-planet model. 
The red dots represent the refined orbital solution with k+1- 
planets as obtained by the recursive periodogram. The resulting 
sampling of the red dots is not uniform in frequency because the 
tested k+1 period is also allowed to adjust. 



While a solution including the 36 and 3.6 days signals provided 
a very extreme reduction of the RMS (from 3.1 to 1.6 m s"' ), 
the peaks in the periodograms of the residuals provided analytic 
FAP estimates too large to be acceptabl e (~ 5%). A qu ick-look 
Bayesian analysis of the same new RVs Tuomil (l2012l) . also in- 
dicated that additional candidates were strongly favoured by the 
data. As we will shown in the analysis section, the RV measure- 
ments of GJ 676A are a textbook example where signal correla- 
tion prevents the detection of lower amplitude signals using pe- 
riodogram methods based on the analysis of the residuals only. 



4. Planetary system : new candidates 

4.1. Recursive periodogram analysis 

For the recursive periodogram analysis, a 1.0 m s"' jitter was 
added in quadrature to the nominal uncertainty of each RV mea- 
surement. This value was chosen because, for any multi-planet 
solution we attempted, about ~ 1.0 ms"' always had to be added 
in quadrature to match the nominal uncertainties to the RMS 
of the residuals. As seen in the top panel of Figure [1] there 



Instead of fitting a trend, we performed a recursive peri- 
odogram search for a second planet with periods between 1.1 
and 50 000 days, obtaining a preferred solution around 4000 
days. The presence of a peak at only twice the time baseline 
indicates the presence of significant curvature (see top-left panel 
in|2]i. As for GJ 676Ab, there is little doubt on the statistical sig- 
nificance of this signal/trend (analytic FAP threshold of 1% is 
around 15, while the signal has an F-ratio of several hundreds). 
As shown in the second panel of Figure [1] the recursive peri- 
odogram (red dots) compared to the periodogram on the residu- 
als (black line) is able to massively improve the significance of 
this second signal thanks to the simultaneous adjustment of the 
orbit of the first candidate. As discussed in the Bayesian analysis 
section, the period and parameters of this candidate are not well 
constrained and only some nominal values are given for refer- 
ence. For detection purposes only, we conservately assume that 
it can be well reproduced by a full Keplerian solution and added 
it to the model. 

After the first two signals were included, the recursive peri- 
odogram search for a third companion revealed one additional 
periodicity at ~ 35.5 days (F-ratio~ 17.5). The analytic FAP was 
0. 156 %, which warranted the empirical FAP computation. In the 
first 1000 trials, 5 trials generated false alarms (FAP ~ 0.5%) 
meaning that more trials are necessary to securely asses if the 
FAP is < 1%. An extended run with 10'* trials produced an em- 
pirical FAP of ~ 0.44% so the candidate was finally accepted. 
The candidate (GJ 676Ae) corresponded to a super-Earth/sub- 
Neptune mass candidate with M sin ; ~ 1 IMe. Even though the 
preferred eccentricity was rather high (~ 0.6) quasi-circular or- 
bits are still allowed by the data. This candidate would receive 
~2.6 times more stellar radiation than the Earth receives form 
the Sun so it would have a hard time keeping oceans of liquid 
water on its surface (ISelsis et alJl2007b . 

Again, this 35.5 days candidate was included in the mod- 
els as a third full Keplerian signal and a recursive periodogram 
search was obtained to look for further companions. A strong 
isolated peak (F-ratio 19.5, P=3.6 days), was the next promising 
one signal, showing an analytic FAP as low as 0.015 %. Only 
1 test over the first 1000 trials generated an spurious peak with 
higher power, indicating that the FAP is significantly lower than 
1% and the candidate was accepted right away. The new can- 
didate (GJ 676Ad) has a minimum mass of ~4.5 and it is 
certainly too close to the star to support liquid water on its sur- 
face. 

The recursive periodogram search for a fifth signal show that 
the next tentative periodicities have analytic FAP at the 10% or 
higher level, which did not satisfy our preliminary detection cri- 
teria so we stopped searching for additional candidates. Even 
though 4 planet signals might seem a lot given that only 75 RVs 
were used, let us note the the amplitudes of the close-in low mass 
companions are relatively high (2 - 3 m s"' , see Figure|2]l com- 
pared to the final RMS of the solution (1.6 m s ' ) and the nom- 
inal uncertainties. 



4.2. Bayesian analysis 

As already discussed, there no doubt that RV data of GJ 676A 
contains the signal of a massive planet (nip - 4.9Mjup) with an 
orbital period of roughly 1060 days (Forveille et al. 201 1) and a 
long period trend. A model with a Keplerian signal and a linear 
trend was chosen as the starting point of the Bayesian analyses. 
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Table 2. Relative posterior probabilities of models with k — 
1,...,4 Keplerian signals {Mk) with or withour a linear trend 
(LT), the Bayesian evidences P{d\Mk), and RMS values. 



k 


P{Mk\d) 


log P{d\Mt) 


RMS [ms-'] 


1+LT 


1.0x10-" 


-224.0 


4.59 


2 


1.8x10^" 


-192.8 


3.00 


3 


3.6xlO-'2 


-179.9 


2.20 


4 


~ 1 


-152.9 


1.67 



While spotting the signature of the massive planet in the RVs 
was trivial, we observed that instead of a linear trend, the sam- 
plings preferred a second Keplerian instead, indicating the pres- 
ence of significant curvature. Therefore, the second model to be 
tested contains the trend modeled as a Keplerian. However, be- 
cause the long-period signal could not be constrained, we fixed 
its eccentricity and period to their most probable values in the pa- 
rameter space (for period, this space was the interval between 1 
and \QTobs days, where Tohs is the baseline of the data) through- 
out the analyses. Because the orbit is only partially covered by 
the time-baseline of the observations, we could not constrain its 
other parameters much either so only the MAP values for this 
candidate are given in Table |3] as a reference. Fixing period and 
eccentricity, however, allowed us to draw representative samples 
from the parameter space and to calculate reliable estimates for 
the Bayesian evidences of each model. The curvature in the long- 
period trend was so clearly present in the data that including cur- 
vature (through a fixed period-eccentricity Keplerian) increased 
the model probability by a factor of l.SxlO''' and decreased the 
RMS of the residuals from 4.59 to 3.00 ms ' (Table|2]i. 

We continued by adding a third Keplerian signal to the sta- 
tistical model and performed samplings of the corresponding pa- 
rameter space. The Markov chains quickly converged to a solu- 
tion that contained the same periodic signal at 35.4 days also 
spotted by the recursive periodograms. The model with ^ = 3 
Keplerians was found to have a posterior probability 2.0x10^ 
times that of a model with k - 2 Keplerians (Table|2]i. The signal 
at 35.4 days corresponds to a planet candidate with a minimum 
mass of 1 1.5 Mg. When sampling the parameter space of a four- 
Keplerian model, we identified a fourth strong signal in the data 
with a period of 3.60 days. This was, again, the same 4th period 
spotted by the recursive periodogram. Our solution of the model 



with k - A further increased the model probability by a factor of 
2.8x10" compared to a model with k - 3, so we could conclude 
that this 3.60 day periodicity was also very confidently present 
in the data (Table |2]i. 

The search for additional periodic signals failed to identify 
significant periodicities so we conclude that the model proba- 
bilities imply the existence of 4 Keplerian signals : the massive 
companion GJ 676Ab at 1.8 AU; a trend with some curvature 
suggesting the presence of another massive giant planet in an 
long period orbit; and two previously unknown planet candidates 
with orbital periods of 3.60 and 35.4 and minimum masses of 
4.4 and 11.5 M® (Tabled Fig.|2|. These signals satisfied all the 
aforementioned detection criteria as the radial velocity ampli- 
tudes were clearly strictly positive and their periods, apart from 
the long-period signal, were well constrained and had narrow 
distributions in the parameter space. In addition to the MAP pa- 
rameter estimates, standard errors, and 99% BCSs in Table[3] we 
show the distributions of the periods, RV amplitudes, and eccen- 
tricities in Fig. [3] These distributions show that - apart from the 
eccentricities of the two new low-mass companions, that peaked 
close to zero - all densities were close to Gaussian shape. 

4.3. Robustness of the Bayesian solution 

To assess the reliability of our solution to the GJ 676A RVs, and 
indeed that of the Bayesian methods in general in assessing the 
existence of Keplerian signals in RV data, we performed a test 
analysis of the first 50 epochs only. The purpose of this test was 
to investigate whether we could spot the same signals, and re- 
ceive the same solution from a smaller number of observations. 
The 50 first epochs have a baseline of approximately 11 99 days 
(roughly two thirds of the full baseline of 1794 days), and be- 
cause of their lower number, we expected them to constrain the 
model parameters less, i.e. yielding broader posterior densities, 
and that the model probabilities are less strongly in favour of 
- possibly even against - the existence of the two new planet 
candidates reported in this work. 

Again, we started with a model containing a single Keplerian 
signal and a linear trend. These were easy to spot from the par- 
tial RV set and we could identify the same massive planet can- 
didate and trend reported by Forveille et al. (2011, , 69 CCF 
measurements were used in that work). However, when we sam- 
pled the parameter space of a two-Keplerian model, we rapidly 
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Fig. 3. Distributions estimating the posterior densities of orbital periods {Px), radial velocity amplitudes {Kx) and eccentricities (e ^ ), 
and three constrained Keplerian signals. The solid curve is a Gaussian density with the same mean (//) and variance (cr^) as the 
parameter distribution has. Additional statistics, mode, skewness {^^) and kurtosis (/i^) of the distributions are also shown. 



Table 3. Orbital solution of the three innermost companions of HD GJ 676A and the excess RV jitter MAP estimates, the standard 
errors, and the 99% BCSs. 



Parameter d e b c(trend)* 

P [days] 3.6000±0.0008 [3.5978, 3.6022] 35.37±0.07 [35.10, 35.45] 1050.3±1.2 [1046.9, 1053.7] 4400~ 

e 0.15±0.09 [0,0.42] 0.24+0.12 [0, 0.56] 0.328+0.004 [0.318, 0.339] 0.2 

/^[ms-'] 2.30+0.32 [1.35, 3.19] 2.62+0.32 [1.66, 3.57] 117.42+0.42 [116.18, 118.66] 41 

w[rad] 5.5+1.9 [0,2;r] 5.8+2.2 [0, 27r] 1.525+0.012 [1.491, 1.557] 6.21 

Mo[rad] 4.1 + 1.7 [0, 2;r] 0.9+2.0 [0, 2;r] 0.957+0.036 [0.844, 1.056] 3.1 

aj [ms-'] 1.38+0.18 [0.95, 1.97] 



Derived parameters 



a [AU] 




0.0413+0.0014 [0.037, 0.045] 


0.187+0.007 [0.17, 0.21] 


1.80+0.07 [1.62, 1.99] 


5.2 


nip sin i 


[Me] 


4.4+0.7 [2.4, 6.4] 


11.5+1.5 [6.5, 15.1] 


1570+100 [1190, 1770] 


951 


wip sin ; 




0.014+0.002 


0.036+0.005 


4.95+0.31 


3.0 


S/S^ 




48.1 


2.3 


0.025 


0.003 



Notes. *** Since all parameters are poorly constrained, the the MAP solution is provided for orientative purposes only. 
Stellar irradiance S at the planet's orbit divided by the flux received by the Earth from the Sun (So). 



discovered another Keplerian signal at 35.5 day period. The 
corresponding two-Keplerian solution together with the linear 
trend increased the posterior probability of the model by a factor 
of 1.0x10'*, which clearly exceeded the detection threshold of 
150. Furthermore, we also identified a third periodicity at 3.60 
days when increasing the complexity of the statistical model by 
adding another Keplerian signal to it. This model was 5.0x10^ 
times more probable than the model with k - 2, s,o we could 
conclude that, three planet candidates and a linear trend were 
already strongly suggested by these initial 50 RVs. Moreover, 
the two new low-amplitude periodic signals satisfy our detec- 
tion criteria by having amplitudes strictly above zero (2.27 [ 1 .00, 
3.41] ms-i and 2.83 [1.48, 4.04] ms"' for GJ 676A e and d, re- 
spectively) and well constrained orbital periods (3.6000 [3.5963, 



3.6027] days and 35.48 [35.16, 35.90] days, respectively). This 
solution is consistent with the one received for the full data set 
in Table |3] which implies that the two new planets could already 
have been detected in the HARPS RVs when the 50th spectrum 
was obtained back in October 2009, possibly even earlier. 

We performed the recursive periodogram analysis on the 
same 50 epochs. Again, the massive GJ 676Ab and the trend 
were also trivially. Then we attempted a recursive periodogram 
search for a third Keplerian. This search spotted the 35.5 days 
signal as the next most likely periodicity, but provided an ana- 
lytic PAP of 15 % which did not satisfy our preliminary detection 
criteria (analytic FAP > 1%). In order to check if the 3.6 days 
candidate could be inferred by periodogram methods, we added 
the 35.5 days signal to the model and performed a recursive pe- 
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riodogram search for a fourth candidate. Although a peak at 3.6 
days peak was present, it was not, by far, the most significant 
periodicity suggested by the recursive periodogram. 

This resuh implies that Bayesian methods are clearly more 
sensitive in detecting low-amplitude signals compared to clas- 
sic periodogram approaches (even compared to our newly de- 
veloped recursive periodogram method). Even if a reasearcher 
prefers to obtain frequentist confirmation (e.g., empirical FAP) 
of a signal before announcing it, early Bayesian detections can 
be used to optimize observational strategies and sample the pe- 
riods of interest. We are conducting simulations with synthetic 
dataset to identify failure modes of the proposed Bayesian meth- 
ods (e.g., identify situations that could generate false positives) 
and refine the detection criteria accordingly. 

5. Conclusions 

We re-derived high precision radial velocities on the public 
HARPS observations of GJ 676A using our newly developed 
HARPS-TERRA software, obtaining a significant improvement 
over the RVs obtained using the CCF approach. We developed 
a recursive periodogram method to enhance the sensitivity of 
least-squares solvers to low amplitude signals in the presence 
of strong multiplanet correlations and provide a recipe to de- 
rive empirical FAPs. We compare the results obtained on the RV 
of GJ 676A to the candidates identified by a Bayesian analy- 
sis obtaining compatible detection of the same four signals. We 
provide the favoured 4 planet solution together with the allowed 
parameter intervals as derived from the Bayesian MCMC sam- 
plings. 

While it is clear that Bayesian methods are more general 
and provide a more complete description of the data, frequen- 
tists methods (e.g., empirical FAP computations) allow a sim- 
pler interpretation of the significance of a detection. The com- 
bination of criteria from both approaches provides great confi- 
dence in our results. We have shown that after the early Bayesian 
detection of 4 planets, 25 more measurements we sufficient to 
confirm the same candidates with periodogram based methods. 
Even if a researcher prefers frequentist confirmation of candi- 
dates, early B ayesian detecti on can be used to optimize follow- 
up programs (lGregorvll20()5i ). In overview, we have shown that 
the confluence of recent data analysis developments (HARPS- 
TERRA, Bayesian toolbox, advanced periodograms) achieve a 
significant boost in sensitivity to very low-mass companions in 
existing datasets. Compared to the significant investment done 
in hardware development, development of more optimal data- 
analysis methods comes at a significantly lower cost thus en- 
abling a more efficient utilization of the observational resources. 

Concerning the new two planet candidates, we find that they 
are both in the sub-Neptune mass regime. The shorter period 
candidate (GJ 676Ad) has a significant probability of transit 
(~ 5% according to Charbonneau et al. 2007) thus encourag- 
ing the photometric follow-up of the star The presence of a long 
period companion (massive planet or a brown dwarf) is now ob- 
viously detected through significant curvature in the trend and a 
period of ~ AOOOdays or longer is preliminary suggested by the 
data. After GJ 876(Riveraet al. 2010) and GJ SSUMav oretal] 
l2009l) . GJ 676A becomes the third M dwarf with 4 planet can- 
didates. Except for the Solar system itself, this planetary system 
has the broadest range of minimum masses and periods reported 
so far (from 5 M® to 5 Mj„p, and from 3.6 days to 4000 or more 
days). Despite of the abundance of candidates, the periods (and 
corresponding semi-major axis) are spaced enough that we do 
not anticipate major dynamical stability issues. Compared to the 



more dynamically packed GJ 581 and GJ 876 systems, the or- 
bits of the candidate planets leave ample room to detect further 
potential candidates in intermediate orbits whenever additional 
RV observations become available. Due to the proximity of GJ 
676A to our Sun (~ 16.4 pc), the long period, ma ssive candidates 
are attractive targets for direct imaging attempts jLagrange et al.l 
|201P). Given that the make-up of stars in binary systems should 
be similar, it would be very interesting to investigate whether GJ 
676B (M3.5V) has been as prolific as GJ 676A in forming all 
kinds of planets. 
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Table 4. Differential HARPS-TERRA RV measurements of GJ 
676A measured in the Solar System Barycenter reference frame. 
Secular acceleration has been subtracted to the RVs. 

BJD RV (7-Rv 

(days) (ms"') (ms"') 
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